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Abstract. The paper concerns numerical algorithms for solving the Beltrami equation fz{z) = 
fi{z)fz{z) for a compactly supported fi. 

First, we study an efficient algorithm that has been proposed in the literature, and present its 
rigorous justification. We then propose a different scheme for solving the Beltrami equation which 
has a comparable speed and accuracy, but has the virtue of a greater simplicity of implementation. 
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1. Introduction. The of existence of a quasiconformal mapping g with a given 
complex dilatation is one of the central problems of classical complex analysis. 

We will first recall a definition of a quasiconformal map. Given an open set Q d C, 
a map / : /(^) C C is said to be absolutely continuous on lines (ACL) in a 

rectangle R d with sides parallel to the x and y axes if / is absolutely continuous 
on almost every horizontal and vertical line in R. The map / is ACL on fi, if it is ACL 
on every rectangle in fl. Partial derivatives fz — {fx — ify)/'^ and = {f^ + ify)/2 
of such map exist a.e. in fJ. 

Definition 1.1. A homeomorphism / : — s- /(fi) is K- quasiconformal if and 
only if the following holds: 

i) f is ACL on every rectangle in f2, 
ii) m < f^lM a.e. tnQ. 

The complex dilatation of a quasiconformal homeomorphism / defined a.e. is 

t^fiz) = fz{z)/fz{z). 

A mapping / with a prescribed complex dilatation /x solves the following Beltrami 
equation: 

(1-1) /.W - mW/.(^)- 

The question of existence of a solution of (|1.1() is addressed in the famous Ahlfors- 
Bers-Boyarskii theorem (see and 0): 

Theorem 1.2. (Ahlfors and Bers, Boyarskii). Let fi G Loc{C) satisfy \\fJ-\\^ < 
K < 1. Then there exists a unique map : C — > C such that f^ satisfies the Beltrami 
equation fg{z) = n{z)f^{z), f^ fixes 0, 1 and oo, and is a (1 + i^)/(l — K)- 
quasiconformal map. 

According to the classical proof of this Theorem, a solution of the Beltrami equa- 
tion can be constructed with the help of two operators: the first operator being the 
Hilbert transform 

(1.2) T[/.](z) = f hm / / d^-Ade, 
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the second — the Cauchy transform 

If the differential /i is compactly supported, then the quasiconformal homcomor- 
phism / is constructed in two steps. First, one solves the following equation: 

(1.4) h* ^ T[fih*] + T[fi]. 

If h* satisfying 11.4|l is found, then / solving is given by 

f{z) = P[f,{h* + l)]{z) + z. 



The operator P is well-defined and is Holder continuous for h E LP{C), p > 2. 
This implies that a solution to H1.4|l must also be constructed in Lp{C), p > 2. The 
celebrated Calderon-Zygmund inequality (cf [3) states that 

\\T[h]\\p < Cp\\h\\p 

for all h G Lp{C) where Cp is some constant, depending on p only, which can be 
chosen so that Cp 1 as p — > 2. While the optimal value of Cp for all p is not known, 
it is known that C2 = 1, and T is an isometry in L2{C): 

\\T[h]h^\\hh. 

To solve ()1.4|l one has to choose a p such that KCp < 1. Then the operator 
h I— > T[^h] is a contraction in Lp(C) and iterations 

(1.5) = T[/i/i"] + r[^] 

converge to a solution of H1.4|l : h* as n ^ oo. 

We can now state the extended version of the Ahlfors-Bers-Boyarskii theorem 
for compactly supported differentials: 

Theorem 1.3. Let fi G Loo(C) be compactly supported and satisfy \\^J■\\^ < K < 
1. Then for every p > 2 such that CpK < 1, the operator h — > T[^{h + 1)] is a 
contraction on Lp(C) with a unique fixed point h* . Moreover, the solution of the 
Beltrami equation, f^, is given by 

(1.6) /'^ - P[^i{h* + 1)] + id, 

and is such that /(O) ~ 0, f is continuous, has distributional derivatives, and fz^l& 
Lp{C). 

We will concentrate on the numerical implementation of the procedure described 
in this Theorem. This is an important practical problem which has been addressed 
by various authors. 

Efficient algorithms for computing the two required transforms T and P have 
been proposed in [H| and ^U]- We will proceed with a very brief description of 
these algorithms. 

As we have pointed out, the Cauchy transform of an Lp-function, p > 2, is well- 
defined and is known to be Holder continuous with exponent 1 — 2/p (cf PP, [3]). 
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In |Sj an effective algorithm for the numerical computation of the following integral 
transform was suggested: 

A comparison with immediately shows that, 

(1.8) p[h]{z) = p[h]{z)-p[hm. 

Represent h and P[h] as a Fourier series: 



(1.9) h{re'')^ hk{ry''', 

k— — oo 
oo 

(1.10) P[h]{re^<') = Pkir)e"'', 

k— — oo 

where the coefficients of the P-transform are given by 

(1.11) p,{r) ^l-j\-^'^<^P[h]{re^<^) d0, 



then, according to |H], 



(1.12) Pk{r) 



2^ Q hk+i{p)dp, k<0, 



OO / \ ^ 

r 



-2 / ( - ) hk+i{p)dp, k>0. 



Next, to obtain similar formulae for the Hilbert transform, assume that h is 
a Holder continuous function compactly supported in an open disk around zero of 
radius R, B{0, R) C C. The Hilbert transform of such function is well-defined (see 
next Section). Represent the Hilbert transform as a Fourier series 



(1.13) T[h]{reJ'^) = V Cfe(r)e''=^ c,{r) = / e-''=^T[/j](re'«) d9. 

In [H] and ^H] the authors arrive at the following expressions for these coefhcients: 

(1.14) co(0) = -2 lim / ^^^^dp, and Cfc(O) = 0, whenever k^O, 

(1.15) Ck{r) = Ak J --j^hk+2{p)dp + Bk J --p-^-hk+2{p)dp + hk+2{r), 
where 

nifiA A / 0,A:>0, r-2(fc + l), fc>0, 

^^- = l2(fc + l),fc<0, ^-dS. = | 0^ 

Given values of h, for instance, on a circular N x M grid that contains the com- 
pact support of one can use a fast Fourier transform (FFT, cf 12 ) to find the 
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values of the coefficients h^. at the radu r^, 1 < i < M. Next, one can use these 
values to construct a piecewise constant, a piecewise linear or a spline approximation 
of the functions hk (the choice of approximation, of course, depends on the known or 
expected smoothness of hk)- This allows one to compute integrals in (|1.12|l and in 
(|1.14() - (|1.15() . Armed with these implementations of the Hilbert and Cauchy trans- 
forms, one can try to solve the Beltrami equation 1)1. first by running iterations 
()1.5(l for some time, and, finally, applying (|1.6(l . It is convenient to use the point-wise 
multiplication of grid values of h and fi + 1 inside the Hilbert transform in l|1.5|l . 
rather than the multiplication of their Fourier series: The order of the computational 
complexity of the point- wise multiplication is 0{NM), as opposed to 0{NM^) for 
the series. The transition from the representation of ft. as a Fourier series to point 
values at each iteration step can be performed with the help of the FFT. This way, 
the computational complexity of one iteration step becomes 0(A^M log2 M). 

The objective of our paper will be two-fold. First, we will present a rigorous math- 
ematical justification of the formulae (|1.14|) - (|1.15|) for the singular operator T different 
from that of [lUj . Second, we will propose an algorithm which avoids programming 
formulae 1)1. 14(l - (|l. 15(1 altogether. It should be noted, that in our experience writing 
the code for the Hilbert transform can be rather cumbersome, therefore this second 
algorithm can be viewed as an improvement. 

2. A Justification of P. Daripa and D. Masiiat's algoritiim for tiie 
Hilbert transform. The Hilbert transform is a singular integral transform, there- 
fore, first of all, one has to specify the conditions under which this integral can be 
computed. The issue of existence of a rather general class of singular integral trans- 
forms (including the Hilbert transform) has been addressed in the classical work |Zj 
of A. P. Calderon and A. Zygmund. The sufficient condition for the existence of the 
Hilbert transform itself can be stated in a very concise way (cf. lllp: 

Lemma 2.1. If h is Holder continuous, then the Hilbert integral transform of h 
exists as a Cauchy principle value. 

In this light, the main result of ^U] can be stated as follows: 

Theorem 2.1. IfT[h]{z) exists as a Cauchy principle value, than the Fourier 
coefficients of this Hilbert transform are given by formulae jJ. J5)) . 

The proof of this Theorem presented in piU starts with an observation that ac- 
cording to H1.2|l and 1)1. 13f) . coefficients of the Hilbert transform are given by the 
following formula: 

(2.1) c,(r) = r e-^^' hm / / ^"1,, pdpdadO. 

p. Daripa and D. Mashat proceed in their proof by computing the Hilbert trans- 
form in the problematic region (p close to r) explicitly. A crucial step in the ensuing 
proof is the interchange of the order of integration in H2.1|l . We feel that the possi- 
bility of this interchange has to be carefully explained: Although the argument that 
we present in this Section is not exceptionally difficult, it nevertheless involves some 
reasoning and has to be carefully carried out. 

We will now proceed with our argument. Let h be Holder continuous with ex- 
ponent 7 and constant A, and let it be compactly supported in an open disk around 
zero of radius R, B{Q, R) C C. Notice, that the variable z in 1)1. 2|l can range over all 
of C and not only over B{0, R). We will start by showing that the Hilbert transform 
of a Holder continous function is bounded. 
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Choose a small 6 < 1/R. For all e < 6R define 



^^I Ie 



B{0,R)\B(z,e) ~ 

HO 



B(0,R)\B{z,5R) 



+ 



HO 



d^AdS, 



(2.2) 



27r J J B{z,SR)\B{z,<i) {£. - 



d^Ad^ 



where r = \z\. 

Lemma 2.2. There exists a constant C such that \T^[h]{z)\ < C for all < e < 
SB and all z € B{0,R). 
Proof. 



\T![h]{z)\ = 



I Ibi. 



HO 

271" J J B{z,5R)\B{z,e) " 
SR p2-K 



d^Ad^ 



Hz + pe'') , 

P e 



£ JO 



pe 



2i9 



dpde 



SR 



£ JO 



^ h{z + pe'(^+-)) -h{z + pe*(«+3-/2)) 



dpde 



< 



|/i(z + pe*^) - h{z + pe 



dpde 



1 f^" p \h {z + pe'(«+^)) -h{z + pe'(^+3V2)) I 



£ JO 



< 



1 f^^ Ap^'ll- e' 



dp + 



1 f^^ Ap^f \l - e'i 



I 



dpde 



dp 



1 



< A\l - - {6Ry = c, 
7 



where we have used Holder continuity. □ 

Next, we will derive the expressions for the Fourier coefficients of the Hilbert 
transform separately for z G intB(0, -R)-{0}, z e dB{0, R),z = and z e C\B{0,R). 

2.1. Case z G intB(0,i?) - {0}. 



Corollary 2.2. 
1 



ipe^c. _ ,,iey dedapdp 



27r2 e^oJ^^Jo 



R ,-2 



h{pe' 



ike 



(pe'« - re'^)2 



dedapdp + hh+2{r) 
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for all < r < R. 

Proof. Let z = re^^ . Recall, that T,[h] = T^\h] + T'^[h] where T\h] and T'^[h] 
are as in (|2.2|l . The first of these operator is clearly bounded, the second is bounded 
by the preceding Lemma. Therefore, |re[/i](z)| is bounded by a constant function, 
while the function = e~*^Te[/i](re**) is bounded by const • e~*^ which is cer- 

tainly in Li([0, 27r)). A direct application of the Lebesgue Dominated Convergence 
Theorem shows that the limiting procedure and the integration over the circle can be 
interchanged: 



/.27r ■ ^27r 

I . ^ 



c,M = — ihn^.Wd^.^hmy^ Mo)de 

^-OJO J JB(Q,R)\B(z.e) ' ' ' re''^)^ e^O ' 



We will denote O = B{0, r + e) \ {B{0, r — e) U B{z, e)) for the purpose of nota- 
tional brevity . Then, for an arbitrary but fixed e one can interchange the order of 
integration (according to Fubini's Theorem) in the following way: 

■ I-2TT p p -ike 
47r2 7o J JBiQ,R)\BiQ,r+e) " rC^'V 



A 2 I II MOtt ^ d^Ad^de 

I-2TT I-2TV ^-ike 



47r2 7o J Jn'^^^^W^^r 



(2-3) + 7^ / / / M07n«W '^C A rf^d^?. 



Next, define Q2 ~ B{z, e'^) \ B{z, e) for some real positive /3 < 1 to be specified 
later {B{z,e) C B{z,e^) since e is always less than 1 by our choice of 6), and fii = 
f2 \ f22 • Using Holder continuity once again, one can show that the difference 

can be bounded in the following way: 

<^ff p^-^ dapdp <^ f f dapdp +^ f f p''^^ dapdp 



27r J Jq 2tt J 7oi Stt j jj^^ 



Ae/3(7-2)A,ea(r!i) + -{e^'' - e^) < AS'<-^)47rre + -( 
27r 7 27r 7 



<2.4re/'(^-2)+i + -(e'3^-e'^). 
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Notice, that 



(2.5) 



lim /, = 



as long as (3{2 — 7) < 1. Now, it is not hard to show that 

d^Ade 1 (r-e)2 



4^ 



•-^■^^ 47r2 J (e - z)2 27r z2 

(see 2ni)j and, therefore, combining (|2.5|) and (|2.6|) . that 



— ^ Urn 



2ir 



/i(Oe- 



1 



J Jn 



27r 



271- 



— - Um 

47r2 e^O 



/i(Oe- 



1 



Um 

e— 



lim 



4^ 



J Jn 

2tt 



d£, A dd£,de - — lim / h{z)e''^'^ - — ^ dO 



, (r - e)2 



2tT e-+0 Jg 



1 



A d^de - — / h{z)e 



J Jn 

2tt 



4^^ 7o 

< lim I, = 0, 



/i(Oe- 



2ir 



d6i 



z^ 

H^)£!!l d^Ad^do 



where in the first equality we again used the Lebesgue Dominated Convergence The- 
orem. We conclude that 



(2.7) —r lim 



(^ — zy ZTT 



/i(z)e-*'=% de. 

Z'' 



Observe, that 

<.27r 



(2.8) 



1 
2^ 



1 







hiz)e-'^^^ de 

z^ 27r 



27r „-ik0 V^oo 



o2ie 



^6* = hk+2ir)- 



The claim of the Corollary follows from (|231), l|2!ZIl and □ 



Now, consider 



2tt ^-ike 



de 



2tx ^-ike 



de 



(^ — re )2 ''^ Jo (^/r — e'^)2 ir'^ Jjiw^^ {S, / r — w) 

We will first compute this integral in the case when |^|/r > 1: 



dw. 



I{i,r) = ^Res 



ui=0 



where 



Bk — 



-2(fc + 1), A: > 0, 
0, fc < 0. 



Similarly, for |^|/r < 1, 
1 



/(C,r) = —Res 



w'^+i(^/r — u')2 



2-K 
H — ^Res 



tu=0 '> 
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where 

Finally, we obtain: 



, fc > 0, 
2{k + l) ,k< 0. 



CfcW = 2:^y^ Kpe"")-^ pdpda + — J hipe'")-^ pdpda + hk+2{r) 



27r Jo Jo pk+ie^(k+2)a 27r 7^ Jo pfc+ie»(fc+2)a 

+ ^A;+2 (r) = Afc / —j—^hk+2{p)dp + Bk I —j—^hk+2{p)dp + hk+2{r). 

Jo P Jr P 

2.2. Case z E dB{0,R). We will briefly outline a counterpart of Corollary 12.21 
for \z\ = R. 

Corollary 2.3. 

1 pR-e p2w p2-n -ike 

-1-1. / / ■, / I t- 



c.(i^) ^-^linry^ X '^'^^"Vo ipe^^ - Re^^Y '''^''"P- 

Proof. According to Lemma 12.21 we can exchange the order of the integration 
and the limit: 



^k{R) = / lini / / H^) di A d^de 



-lim / / / h{0 ,, d^Ad^de = \imcUR). 

An^e^Oj^ J JB(0,R)\B(z.e) {£. - Re'^Y 



Denote n = 5(0, R) \ (5(0, i? - e) U B{z, e)). Then, we obtain for an arbitrary 
but fixed e: 



R-£ i-2-K 



= Jo '^''^"^Jo ipe--Re-Y 

■ p2t! p p -ike 

^^L J 

Again, define il2 — B{z,€^) \ B(z,e), for some real positive /? < 1, and fii = 
\ . As in Corollary 12.21 one shows that the difference H2.4|l can be bounded by 

^^g/3(7-2)+i (^/27)(e^7 _ £7)^ and that lim.^o-^e = as long as (3{2 - 7) < 1. 
However, h{z) = for all z E dB{0, R), hence 



and therefore 



I 11 dihdide 



4^' Jo J Jn (e - z) 



(2.10) ^ li- y '^^^w^r;^ ^ - 0- 

The claim of the Lemma follows from (|2.9|l and H2.10|) . □ 
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To obtain an expression for Ck{R), consider the integral J^^ er"^^^ j — i?e*^)^ dB. 
It is not hard to show that this integral is equal to ~ Akir R'' / S,''^'^ , and, therefore, 

c.(i?) = -^limy^ X ipe^o. _ R,.ey dOdapdp 



. . h{pe''')Ak-j—!j dapdp 

^7'' Jo Jo ? " 



= ^fey -^^hk+2{p) dp. 

2.3. Case z = 0. We will prove a counterpart of Corollary (|2.2|l at the same 
time obtaining an expression for Cfc(O). 
Corollary 2.4. 

f 0, fc^O, 

Proof. As before, according to Lemma l2.2l for an arbitrary but fixed e one can use 
the Lebesgue Dominated Convergence Theorem to change the order of the integration 
and the limit: 

47r^ 7o <'^°J JB(0,R)\B{0,e) K 



hm / / / d^AdCde 

47r2 £^o7o J JB{OM)\B{0,e) C 



1 f" /•'"E:=-oo/^n(p)e™" 



27r2 Jo Jo P^e2«" 
0, A: 7^ 0, 

io p 



dOdapdp 



-2C!^dp, k = 0. 



The Holder property of h and the fact that ft.2(0) = show that the last integral 
(fc = 0) is convergent. □ 



2.4. Case z G C \ B{0, R). In this case the limiting procedure and the angular 
integration in 12.111 can be interchanged right away. 



^ de = ^Res 



/o (C-re^^) 
and 



w=0 S, 
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27r 

A,., 



B(Q,R) 



inO 



-/c+2 



pdpdO 



n— — oo 



p- 



,fc+i 



i'k+2 



(p) dp. 



3. Recursive formulae. Define 

dk{r,r')^2{k + l) f -(-] hk+2{p)dp, 

Jr P\P J 

then, using (|1.14|l - p.l6|l . one readily gets the following recursive formula for the 
coefficients of the Hilbert transform: 

(3.1) Cfc(r) - (^)' (cfe(r') - hk+2{r') - dk{r, r')) + hk+2{r), 

which holds for any pair of nonnegative radii r and r' for which the ratio {r/r')^ is 
defined. 

We will also note, that it is convenient to compute Pfe(r) in (|1.12|) using the 
following recursive relations: 

(3.2) pk{r)^(^^}j%k{r')~ek{r,r'), ek{r,r') = 2 {^-^ hk+i{p)dp. 

These recursive formulae are extremely useful because they allow to compute the 
value of a coefficient at the i-th radius using its value at either {i + l)-th, or [i — l)-th 
radius, thus avoiding the full integration in formulae (|1.12() and H1.14|l - H1.16|l . 

For reference purposes we will rewrite formulae 1)3. l|l and (|3.2I) for a piecewise- 
constant extrapolation of data hk(ri) in Appendix A. Strictly speaking, our derivation 
of the expressions for the Fourier coefficients of the Hilbert transform does not hold for 
a discontinuous h with piecewise-constant coefficients. However, given an e > there 
always is a Holder-continous function whose coefficients are e-close to hk (r) in the Li 
norm, and, therefore, there always is a Holder-continous smoothening of piecewise- 
constant data such that the coefficients of its Hilbert transform are arbitrarily close 
to those given by expressions (|5.1|) ^ (|6.3|) in the sup norm. 

4. Another way to compute T. In this Section we propose an alternative way 
for computing T[h]. By a theorem, presented in the following identity holds for 
all h £ Cq (twice differentiable h with a compact support): 

(4.1) T[h]{z) = d,P[h]{z) = d,P[h]{z). 

As we have already mentioned, in 9 Daripa has suggested an effective algorithm 
for a numerical computation of integral transform (|1.7(l . The algorithm in 9 for 
computing P is defined in terms of a circular N x M grid and its computational 
complexity is 0{NM\og2M). In order to use this algorithm it is necessary to find 
an approximation for /i^, given the values of /i on a circular grid. A formula for 
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in polar coordinates is supplied in the following Lemma, whose easy prove will be 
omitted. 

Lemma 4.1. Suppose that f is C^(M?), and in polar coordinates f is represented 
by g, i.e. f{x,y) = g{r,9), where x + iy = re^^ . Then 



(4.2) Mx,y) = —{g,{r,e)+ge{r,0)/ir), 

ie'" 

(4.3) Mx,y) = —{-ig,(r,e)+ge{r,9)/r). 



5. Comparison of the two algorithms. We have implemented the algorithm 
for the Cauchy transform (see formulae H1.7|l - (|1.12|1 and H3.2() ), both algorithms for 
the Hilbert transform f lfni|l - fLT(^ . l|TT|l and (glj) for the piecewise linear approx- 
imation of Fourier coefficients as a set of routines in the programming language Ada 
95 (cf J3] for the language standard). We have used the public version 3.15p of the 
GNAT compiler Our programs can be found at |15|. 

Of particular interest to us will be a comparison of the two implementations of 
the iteration scheme h i-^ T[/i(/i+ 1)]: in one of them, the Hilbert transform has been 
computed through formula 1)3. in the other, as the z-derivative of transform 1)1. 7|l . 
Below, we will refer to these two iteration schemes as Scheme 1 and Scheme 2. 

The prototypical Beltrami differential for algorithm testing in the literature seems 
to be /i(z) — z^z'' in some disk and /i(z) = outside of that disk. Both Hilbert and 
Cauchy transforms of this differential, as well as the solution of the Beltrami equation 
itself, can be found exactly. Here we will quote the formula for the Hilbert transform: 

^ pg{R~\z\)^^,,_^^, ^ 
q + l 

(5.1) + [r,{q + l-p)-r,{R- \z\)] P^±_1rH,+i) .P-.-^ 

q+l 

where rj is the Heaviside step function. The derivation of this formula is relatively 
straightforward, and can be done either by a direct calculation or by differentiat- 
ing the expression for the Cauchy transform of z^z' which have appeared in several 
publications (cf e.g. 8 ). 

This formula suggests yet one more algorithm for the Hilbert transform of Cg 
functions representable as finite power series on its compact support. Let h{z,z) = 
rj{R — IzDY^I^j^Qttijz'^z^ be such function. Then one can first use the linearity of 
the Hilbert transform to represent T[h] as T[h]{z,z) = X]"j"o ^'.i-^I'^*^"']' compute 
each r[z*z^] in this series using (|5.1|l and sum up. We will refer to this procedure as 
Scheme 3. Of course, this Scheme has a limited value, as many differentials do not 
have sufficient smoothness, or are sometimes given in numerical experiments simply 
as an array of values on a grid. 

Here, we will use a quartic test Beltrami differential which readily yields itself to 
the procedure that we just described: 

^'^'^-\ 0, |z|>V2s, 

where a and s are some real positive parameters. Our /i is given by a real-valued, 
non-negative, radially symmetric Cq function on C with a single maximum at zero. 
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compactly supported in the disk of radius \/2s. We would like to emphasize, that 
iterations i— > T[h"{fj,+ 1)] can be performed exactly using Scheme 3 (that is, the 
only error is the rounding error of our computing device). In particular, there is no 
need for any approximation of the functions being transformed, as it was the case with 
the choice of the spline approximation of the coefficients hk in formulae (|1.14(l - (|1.15|l . 
We will use Scheme 3 as a measuring stick for the accuracy of Schemes 1 and 2. 

Let g{r,6) = X^fe 5fc('")s''^^ be smooth, then it follows from 14. 2|) that the Fourier 
coefficients of the z-derivative of this function, denoted by ^, are given by the following 
exact expression: 



(5.3) 



g'kir) + kgk{r)/r 



We have used this formula to compute the coefficients of the z-derivative of the 
iteration h i— > dzP[fi{h + 1)] for our choice (j5.2|l of 5fc(r)'s were supphed by right 
two-point approximations through values of gkir) on the radial grid. 

We should also mentioned that we have used the fast Fourier transform (FFT) 
together with the point-wise multiplication of function values on the grid to perform 
the multiplication of fi and h in both schemes. Indeed, the computational cost of the 
FFT is of order 0{NM log2 M), that of the point- wise multiplication of functional 
values is 0{NM); at the same time, the complexity of the multiplication of two Fourier 
series representing fj, and h would be 0{NA'P). 

We subjected the two schemes to several test on an 1.5 GHz Intel Pentium class 
PC. Our first test was that of convergence of the two iteration schemes. Table 1 
compares the convergence rate (measured by the norm \\h" — /i"~^|joo rounded to 
3 significant figures) for both iteration schemes and several values of n. We have 
used zero initial conditions, h'^ = n m both cases. We have fixed the values of the 
parameters in the Beltrami differential to be a = 0.5 and s — 1.0. The dimensions of 
the radial grid were N = 500 (i.e. the disk of radius \/2 contained 500 radii of the 
grid) and Af = 512. 



number of iterations, n 


5 


10 


20 


50 


Scheme 1 


1.17 X IQ-^ 


3.83 X 10-" 


1.80 X 10-14 


1.17 X 10-19 


Scheme 2 


1.20 X 10"^ 


3.91 X 10-1° 


6.79 X 10-14 


8.27 X 10-20 



Table 1. Convergence rate (||/i" — /i" ^||oo) for the two iteration schemes. 



It is not surprising that the convergence rates of the schemes are comparable: 
These rates should depend only on the value of the essential supremum of and not 
on the particular realization of the Hilbert transform. 

Next, we have measured the execution time (rounded to nearest seconds) of 10 
iteration steps for several grids. Table 2 displays the results (the values of a, s and 
/i" were as before). 



N X M 


500 X 256 


1000 X 512 


5000 X 1024 


7000 X 2048 


Scheme 1 


8 


31 


313 


965 


Scheme 2 


10 


38 


375 


1166 



Table 2. Execution time for several grids (in sec). 
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As it should be expected, Scheme 2 takes somewhat longer: The computational 
complexity of both integral transforms is the same, 0{NM\og2M), however, the 
iteration h i-^ dzP[fJ.{h + 1)] also involves a differentiation whose computational com- 
plexity is 0{NM). Nevertheless, for the tested range of grids, the speed benefits of 
Scheme 1 probably do not justify an extra effort in programming the Hilbert trans- 
form: As we have already mentioned, in our experience writing and testing such code 
might be rather cumbersome and lengthy, and could significantly increase the project 
development time. 

The loo norm of the difference of the outcomes of Schemes 1 and 3 and Schemes 
2 and 3 after 10 steps (with the same parameters and initial conditions as before) for 
several grids is reported in Table 3. As we have already explained it. Scheme 3 serves 
as a measuring stick for our purposes. 



N xM 


500 X 256 


1000 X 512 


5000 X 1024 


7000 X 2048 


\\hl-hl\\oo 


4.00 X 10-^ 


1.00 X 10"^ 


7.12 X 10-** 


1.89 X 10-^ 


1^2 ~ ^il oo 


8.39 X 10-5 


4.20 X 10-5 


8.39 X IQ-^' 


5.99 X 10-6 



Tabic 3. Comparison of the performance of Schemes 1 and 2 against Scheme 3. 

With regards to the accuracy of Schemes 1 and 2, Table 3 demonstrates that 
Scheme 1 is the scheme of choice. A bigger error in Scheme 2 can be attributed to 
the fact that we are using a numerical approximation for the differentiation in the 
iteration h i— s- dzP[^ji{h + 1)]. This approximation introduces an extra error which is 
absent in Scheme 1. However, if the degree of differentiability of h is known, one can 
use a proper spline approximation of the derivative and thus significantly reduce this 
error. We still feel that in those cases when a high accuracy is not required. Scheme 
2 can be more preferable. 

6. Appendix A. We will now return to formula H3.1(l . Consider a collection of 
real numbers = ri < r2 < . . . < r^^i < = R < r^+i < . . . < rjv that specify a 
radial grid. We will denote the midpoint of the interval [ri,ri+i] by p^; go will be 
by definition. 

If /ife's are compactly supported piecewise- constant functions 

hk.i, r e [ft-i, ft)) ^ <'i < L, 
0, r > QL+i, 



hk{r) 



then the coefficients of the Cauchy and the Hilbert transforms are given by the fol- 
lowing set of expressions. 

Coefficients of the Hilbert transform for piecewise-constant data. 
• r = ri 



(6.1) Cfc(ri) = 

• re (ft-1, ft], l<i< L 



Co(r2)-S2(r2)-2s2(r2)ln^, fc = 0, 



0, 



fcT^O. 



(^) i'^k^^i+i) ~ ^k+2{ri+i)) - dk,i{r) + Sk+2{r), k>l, 
(6.2) Cfc(r) =<(co(r,+i)-S2(r»+i)-2s2(r»+i)ln^-2s2(r01nf + S2(r), k ^ 0, 

(n~r) i'^kiri-i) ~ Sk+2iri-i)) + bka{r) + Sk+2{r), k<0, 
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where 



4,i(r) = 2 
bk,i{r) = 2 



k+1 
k 

k + 1 
k 



Sfe+2(r-i+i) 
Sfe+2(n-i) 



• re L <i<N, 

(6.3) Cfc(r) = 



+ SA;+2(r-i) 1 T 



+ Sk+2{ri) I -r 1 



0, fc^^O, 
-Cfe(rj_i), fc < 0. 



Coefficients of the Cauchy transform for piecewise- constant data. 
• r = ri 

_ / co(r2) - 2si(r2)eo - 2si(r2)(r2 - £ii), A; = 0, 
0, k^O. 



Cfe(ri) 

r e (ft-i,^], 1 < i < i, 



^ci(r,+i) - 2r,S2(r,) In f - 2rs2(r,+i) In 
(t^t) c/c(rj_i) +6fe,i(r), 



fc > 1, 
fc = 1, 

fc <= 0, 



where 



dk,i{r) 



2r 
1-fc 

2r 



fe-i 



1 j,k 1 



Sfe+i(?'i) -j-Y - 1 + Sfc+i(r,+i) - 



1-fc 

• e (ft-i,^], L<i<N, 



sk+iin-i) -H - ihr 1 + sk+i{ri) 1 1 



44^ 



c;s(r) 



0, 



k^O, 



-Cfe(ri_i), fc < 0. 
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